BSLMM GTEx results

"%&%" = function(a,b) paste(a,b,sep="")
library(ggplot2)
library(dplyr)
library(tidyr)
my.dir <- '~/GitHub/GenArch/GenArchPaper/BSLMM/'
my.vol <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/BSLMM_exp/'
h2.TS <- read.table(my.dir %&% 'bslmm_gtex_results/GTEx_Tissue-Specific_local_h2_se_geneinfo.txt',header=TRUE, sep='\t')
h2.TW <- read.table(my.dir %&% 'bslmm_gtex_results/GTEx_Tissue-Wide_local_h2_se_geneinfo.txt',header=TRUE, sep='\t')

Cross-tissue by OTD

arrange by median PVE

hun <- read.table(my.dir %&% 'bslmm_gtex_results/cross-tissue_exp_BSLMM-s100K_iterations_all_chr1-22_2015-07-20.txt',header=T)
data <- hun %>% arrange(pve50) %>% mutate(position=1:length(pve50),`pge025>0.01`=pge025>0.01) 
data <- data[complete.cases(data),]
genenames <- h2.TW %>% select(gene=EnsemblGeneID,genename=AssociatedGeneName)
data <- left_join(data,genenames,by='gene')
#head(data,n=3L)
#tail(data,n=3L)
ggplot(data,aes(x=position,y=pve50,ymin=pve025,ymax=pve975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

ggplot(data,aes(x=position,y=pge50,ymin=pge025,ymax=pge975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

ggplot(data,aes(x=position,y=n_gamma50,ymin=n_gamma025,ymax=n_gamma975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

arrange by median PGE

hun <- read.table(my.dir %&% 'bslmm_gtex_results/cross-tissue_exp_BSLMM-s100K_iterations_all_chr1-22_2015-07-20.txt',header=T)
data <- hun %>% arrange(pge50) %>% mutate(position=1:length(pge50),`pge025>0.01`=pge025>0.01) 
data <- data[complete.cases(data),]
genenames <- h2.TW %>% select(gene=EnsemblGeneID,genename=AssociatedGeneName)
data <- left_join(data,genenames,by='gene')
#head(data,n=3L)
#tail(data,n=3L)
ggplot(data,aes(x=position,y=pve50,ymin=pve025,ymax=pve975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

ggplot(data,aes(x=position,y=pge50,ymin=pge025,ymax=pge975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

ggplot(data,aes(x=position,y=n_gamma50,ymin=n_gamma025,ymax=n_gamma975,col=`pge025>0.01`)) + geom_pointrange(col='gray')+geom_point()

check correlations

ggplot(data,aes(x=pve50,y=pge50,ymin=pge025,ymax=pge975,col=`pge025>0.01`)) + geom_pointrange(col='gray') + geom_point() 

cor.test(hun$pge50,hun$pve50)

## 
##  Pearson's product-moment correlation
## 
## data:  hun$pge50 and hun$pve50
## t = 130.5, df = 16948, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7003744 0.7153942
## sample estimates:
##       cor 
## 0.7079643

ggplot(data,aes(x=log10(n_gamma50),y=pge50,col=`pge025>0.01`)) + geom_point(alpha=0.3)

Compare GCTA marginal h2 to BSLMM median PVE (local SNPs)

gcta <- h2.TW %>% select(gene=EnsemblGeneID,genename=AssociatedGeneName,local.h2=h2.CrossTissue)
all <- inner_join(hun,gcta,by='gene') %>% arrange(desc(pve50))

## Warning in inner_join_impl(x, y, by$x, by$y): joining factors with
## different levels, coercing to character vector

dim(all)

## [1] 10571    21

ggplot(all,aes(x=local.h2,y=pve50))+geom_point(alpha=0.4)+coord_cartesian(xlim=c(0,1),ylim=c(0,1))+xlab(expression("GCTA h"^2))+ylab('BSLMM PVE')+geom_abline(c(0,1),color='red')

## Warning: Removed 1 rows containing missing values (geom_point).

cor.test(all$local.h2,all$pve50)

## 
##  Pearson's product-moment correlation
## 
## data:  all$local.h2 and all$pve50
## t = 88.659, df = 10568, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6420295 0.6638981
## sample estimates:
##       cor 
## 0.6530999

check tissue-specific correlations PGE vs PVE

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
ts <- data.frame()
for(tis in tislist){
  data <- read.table(my.vol %&% tis %&% '_TS_exp_BSLMM-s100K_iterations_all_chr1-22_2015-08-06.txt',header=T,sep="\t")  
  subdata <- select(data,gene,pve50,pge50,pge025,pge975) %>% mutate(tissue=tis,`pge025>0.01`=pge025>0.01)
  res<-cor.test(subdata$pge50,subdata$pve50)
  cat(tis,"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  ts <- rbind(ts,subdata)
}

## Adipose-Subcutaneous     Pearson R= 0.441    P-value= 0 
## Artery-Tibial    Pearson R= 0.473    P-value= 0 
## Heart-LeftVentricle  Pearson R= 0.443    P-value= 0 
## Lung     Pearson R= 0.426    P-value= 0 
## Muscle-Skeletal  Pearson R= 0.537    P-value= 0 
## Nerve-Tibial     Pearson R= 0.518    P-value= 0 
## Skin-SunExposed(Lowerleg)    Pearson R= 0.501    P-value= 0 
## Thyroid  Pearson R= 0.572    P-value= 0 
## WholeBlood   Pearson R= 0.558    P-value= 0

##exclude gene with pve ~0 and pge =1
#ts<-filter(ts,gene!="ENSG00000196409.7")
ggplot(ts,aes(x=pve50,y=pge50,ymin=pge025,ymax=pge975,color=`pge025>0.01`) ) + facet_wrap(~tissue,ncol=3) + geom_pointrange(col='gray')+geom_point()+theme_bw()+ggtitle("BSLMM Tissue-Specific")

check tissue-wide correlations PGE v PVE

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
tw <- data.frame()
for(tis in tislist){
  data <- read.table(my.vol %&% tis %&% '_TW_exp_BSLMM-s100K_iterations_all_chr1-22_2015-10-18.txt',header=T,sep="\t")  
  subdata <- select(data,gene,pve50,pge50,pge025,pge975) %>% mutate(tissue=tis,`pge025>0.01`=pge025>0.01)
  res<-cor.test(subdata$pge50,subdata$pve50)
  cat(tis,"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  tw <- rbind(tw,subdata)
}

## Adipose-Subcutaneous     Pearson R= 0.653    P-value= 0 
## Artery-Tibial    Pearson R= 0.66     P-value= 0 
## Heart-LeftVentricle  Pearson R= 0.589    P-value= 0 
## Lung     Pearson R= 0.611    P-value= 0 
## Muscle-Skeletal  Pearson R= 0.662    P-value= 0 
## Nerve-Tibial     Pearson R= 0.671    P-value= 0 
## Skin-SunExposed(Lowerleg)    Pearson R= 0.654    P-value= 0 
## Thyroid  Pearson R= 0.67     P-value= 0 
## WholeBlood   Pearson R= 0.663    P-value= 0

ggplot(tw,aes(x=pve50,y=pge50,ymin=pge025,ymax=pge975,color=`pge025>0.01`) ) + facet_wrap(~tissue,ncol=3) + geom_pointrange(col='gray')+geom_point()+theme_bw()+ggtitle("BSLMM Tissue-Wide")

plot tissue-specific PVE vs GCTA marginal h2

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
tislist2 <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.spaces.tissue.list',sep="\n",what="character")
ts <- data.frame()
for(i in 1:length(tislist)){
  bs <- read.table(my.vol %&% tislist[i] %&% '_TS_exp_BSLMM-s100K_iterations_all_chr1-22_2015-08-06.txt',header=T,sep="\t") %>% select(gene,pve50)
  h2 <- read.table("/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/gtex-h2-estimates/GTEx.resid.tissue-specific.h2_" %&% tislist2[i] %&% "_marginal.local_2015-03-23.txt",header=T, sep="\t") %>% select(tissue,ensid,h2) %>% mutate(gene=ensid)
  subdata <- inner_join(h2,bs,by="gene")
  res<-cor.test(subdata$pve50,subdata$h2)
  cat(tislist2[i],"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  ts <- rbind(ts,subdata)
}

## Adipose - Subcutaneous   Pearson R= 0.565    P-value= 0 
## Artery - Tibial  Pearson R= 0.611    P-value= 0 
## Heart - Left Ventricle   Pearson R= 0.589    P-value= 0 
## Lung     Pearson R= 0.578    P-value= 0 
## Muscle - Skeletal    Pearson R= 0.629    P-value= 0 
## Nerve - Tibial   Pearson R= 0.587    P-value= 0 
## Skin - Sun Exposed (Lower leg)   Pearson R= 0.648    P-value= 0 
## Thyroid  Pearson R= 0.608    P-value= 0 
## Whole Blood  Pearson R= 0.688    P-value= 0

ggplot(ts,aes(x=h2,y=pve50))+geom_point(alpha=0.4)+coord_cartesian(xlim=c(0,1),ylim=c(0,1))+xlab(expression("GCTA h"^2))+ylab('BSLMM PVE')+geom_abline(c(0,1),color='red') + facet_wrap(~tissue,ncol=3)+theme_bw()+ggtitle("Tissue-Specific")

plot tissue-wide PVE vs GCTA marginal h2

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
tw <- data.frame()
for(i in 1:length(tislist)){
  bs <- read.table(my.vol %&% tislist[i] %&% '_TW_exp_BSLMM-s100K_iterations_all_chr1-22_2015-10-18.txt',header=T,sep="\t") %>% select(gene,pve50)
  h2 <- read.table("/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/gtex-h2-estimates/GTEx.tissue-wide.h2_" %&% tislist[i] %&% "_marginal.local_2015-03-24.txt",header=T, sep="\t") %>% select(tissue,ensid,h2) %>% mutate(gene=ensid)
  subdata <- inner_join(h2,bs,by="gene")
  res<-cor.test(subdata$pve50,subdata$h2)
  cat(tislist[i],"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  tw <- rbind(tw,subdata)
}

## Adipose-Subcutaneous     Pearson R= 0.414    P-value= 0 
## Artery-Tibial    Pearson R= 0.414    P-value= 0 
## Heart-LeftVentricle  Pearson R= 0.307    P-value= 0 
## Lung     Pearson R= 0.377    P-value= 0 
## Muscle-Skeletal  Pearson R= 0.459    P-value= 0 
## Nerve-Tibial     Pearson R= 0.442    P-value= 0 
## Skin-SunExposed(Lowerleg)    Pearson R= 0.459    P-value= 0 
## Thyroid  Pearson R= 0.465    P-value= 0 
## WholeBlood   Pearson R= 0.38     P-value= 0

ggplot(tw,aes(x=h2,y=pve50))+geom_point(alpha=0.4)+coord_cartesian(xlim=c(0,1),ylim=c(0,1))+xlab(expression("GCTA h"^2))+ylab('BSLMM PVE')+geom_abline(c(0,1),color='red') + facet_wrap(~tissue,ncol=3)+theme_bw()+ggtitle("Tissue-Wide")